Assignment 4-Detroit

options(java.parameters = "-Xmx6G")

library(r5r)

## Loading required namespace: sf

## Loading required namespace: rJava

## Please make sure you have already allocated some memory to Java by running:
##   options(java.parameters = '-Xmx2G').
## Currently, Java memory is set to -Xmx6G

library(osmextract)

## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
## Check the package website, https://docs.ropensci.org/osmextract/, for more details.

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.4     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(sf)

## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

library(ggthemes)
library(ggspatial)
library(tigris)

## To enable 
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.

library(wesanderson)
library(tidytransit)
library(stars)

## Loading required package: abind

library(magick)

## Linking to ImageMagick 6.9.12.3
## Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fontconfig, x11

dir.create("networks")

download.file("https://www.detroitmi.gov/Portals/0/docs/deptoftransportation/pdfs/ddot_gtfs.zip", file.path("networks","ddot_gtfs.zip"), mode = "wb", quiet=TRUE)

Detroit_file <- oe_match("Detroit")

Detroit_streets <- oe_read("networks/bbbike_detroit.pbf",
                   download_directory = "networks",
                   layer = "lines",
                   quiet = TRUE) %>%
 
  filter(!is.na(highway))

michigan_state_plane <- "+proj=lcc +lat_1=43.66666666666666 +lat_2=42.1 +lat_0=41.5 +lon_0=-84.36666666666666 +x_0=3999999.999984 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048 +no_defs"

Detroit_city_limits <- places("Michigan") %>%
  filter(NAME == "Detroit") %>%
  st_transform(crs = st_crs(Detroit_streets))


Detroit_streets <- Detroit_streets[Detroit_city_limits,]

ggplot() +
  geom_sf(data = Detroit_streets, alpha = 0.1) +
  coord_sf(crs = michigan_state_plane)

rJava::.jgc(R.gc = TRUE)

grid <- st_sf(st_make_grid(Detroit_city_limits, 
                           square = FALSE, 
                           n = c(100,100),
                           what = "polygons"))%>%
  st_filter(Detroit_city_limits) 

colnames(grid) <- "geometry"
st_geometry(grid) <- "geometry"

grid <- grid %>%
  mutate(id = seq(1, length(grid$geometry), by=1))


# ggplot() +
#   geom_sf(data = grid) +
#   theme_map()

#Detroit School Work

Detroit_school <- oe_read("networks/bbbike_detroit.pbf",
                   provider = "openstreetmap_fr",
                   download_directory = "networks",
                   layer = "points",
                   quiet = TRUE) %>%
  filter(str_detect(other_tags,'"amenity"=>"school"')) %>%
  st_filter(Detroit_city_limits) %>%
  rename(id = osm_id)

#Aeshna’s Schools Map, 0A

# ggplot(Detroit_streets) +
#   geom_sf(color = 'azure2', alpha = 0.2) +
#   geom_sf(data = Detroit_school, color = "pink4") +
#   coord_sf(crs = michigan_state_plane)  +
#   theme_void()

plot(image_read("images/Aeshna_0A.png"))

grid_points <- st_centroid(grid)

ggplot() +
  geom_sf(data = grid_points, size = 0.75) +
  geom_sf(data = Detroit_school, color = "pink4") +
  theme_map()

r5r_core <- setup_r5("networks", verbose = FALSE)

ttm_school <- travel_time_matrix(r5r_core = r5r_core,
                          origins = st_transform(x = Detroit_school, crs = "WGS84"),
                          destinations = st_transform(x = grid_points, crs = "WGS84"),
                          mode = "BICYCLE",
                          departure_datetime = as.POSIXct("15-11-2021 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S"),
                          max_walk_dist = 1000,
                          max_trip_duration = 480,
                          verbose = FALSE)

tt_wide_school <- ttm_school %>%
  pivot_wider(names_from = fromId, 
              names_prefix = "from", values_from = travel_time) %>%
  rename(id = toId) %>% 
  merge(grid) %>%
  replace(is.na(.), 999) %>%
  rowwise() %>%
  mutate(from_any = min(c_across(starts_with("from")), na.rm = TRUE))

st_geometry(tt_wide_school) <- "geometry"

#Aeshna’s Cycling to School Map, 2A

# ggplot() +
#   geom_sf(data = tt_wide_school, 
#           aes(fill = from_any), 
#           color = NA) +
#   geom_sf(data = Detroit_school, color = 'ivory3') +
#   geom_sf(data = Detroit_streets , alpha = 0.05)+
#   scale_fill_gradientn(colors = wes_palette(name = "Chevalier1", n= 3, type = "continuous"),
#         name = "Cycling time to\nthe nearest school\n(minutes)",
#         position = "right") +
#   coord_sf(crs = michigan_state_plane) +
#   theme_map()+
#   theme(legend.position = "right")

plot(image_read("images/Aeshna_2A.png"))

#Aeshna’s School Isochrone Map, 1A

# iso_pallete <- wes_palette("Royal2", n = 5)
# 
# iso5min_school <- tt_wide_school[tt_wide_school$from_any < 6,] %>%
#   st_union()
# 
# iso10min_school <- tt_wide_school[tt_wide_school$from_any < 11,] %>%
#   st_union()
# 
# iso15min_school <- tt_wide_school[tt_wide_school$from_any < 16,] %>%
#   st_union()
# 
# ggplot(Detroit_streets) +
#   geom_sf(data = iso15min_school, 
#           aes(fill = "Area within 15 minutes"), 
#           color = NA) +
#   geom_sf(data = iso10min_school, 
#           aes(fill = "Area within 10 minutes"), 
#           color = NA) +
#   geom_sf(data = iso5min_school, 
#           aes(fill = "Area within 5 minutes"), 
#           color = NA) +
#   geom_sf(alpha = 0.1) +
#   scale_fill_manual(values = c(iso_pallete[1], 
#                                iso_pallete[3],
#                                iso_pallete[5]),
#         name = "Cycling \ntime to the\nnearest school\n(minutes)") +
#   coord_sf(crs = michigan_state_plane) +
#   theme_map()

plot(image_read("images/Aeshna_1A.png"))

school_grid <- grid %>%
  mutate(num_school = lengths(st_covers(grid, Detroit_school)))

school_points <- st_centroid(school_grid)

ggplot(school_points) +
  geom_sf(aes(color = as.character(num_school))) +
  scale_color_manual(values = c("lightgray", "skyblue", "deepskyblue4", "blue4", "black"), 
                    name = "Number of\nschools") +
  theme_void()

Detroit_access_school <- accessibility(r5r_core, 
                                        origins = school_points,
                                        destinations = school_points,
                                        mode = "BICYCLE",
                                        opportunities_colname = "num_school",
                                        decay_function = "step",
                                        cutoffs = 11,
                                        max_walk_dist = 5000,
                                        time_window = 120,
                                        percentile = 50,
                                        verbose = FALSE) %>%
mutate(id = as.numeric(from_id)) %>%
  merge(grid)

#Aeshna’s Schools Reachable by Cycling Map, 3A

# st_geometry(Detroit_access_school) <- "geometry"
# 
# ggplot(Detroit_access_school) +
#   geom_sf(aes(fill = accessibility), color = NA) +
#   scale_fill_viridis_c(name = "Schools within \n a 20-minute \nbicycle ride") +
#   coord_sf(crs = michigan_state_plane) +
#   theme_void()

plot(image_read("images/Aeshna_3A.png"))

Detroit_school_access2 <- accessibility(r5r_core,
                        origins = school_points,
                        destinations = school_points,
                        mode = "BICYCLE",
                        opportunities_colname = "num_school",
                        decay_function = "exponential",
                        cutoffs = 11,
                        departure_datetime = as.POSIXct("15-11-2021 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S"),
                        max_walk_dist = 5000,
                        time_window = 120,
                        percentiles = 50,
                        verbose = FALSE) %>%
  mutate(id = as.numeric(from_id)) %>%
  merge(grid)

st_geometry(Detroit_school_access2) <- "geometry"

ggplot(Detroit_school_access2) +
  geom_sf(aes(fill = accessibility), color = NA) +
  scale_fill_viridis_c(name = "Accessiblity score") +
  coord_sf(crs = michigan_state_plane) +
  theme_void()

rJava::.jgc(R.gc = TRUE)

st_write(Detroit_school_access2, 'Detroit_school_access.geojson', append=FALSE, quiet=TRUE )

access_poly_school <- st_read("Detroit_school_access.geojson", quiet=TRUE)

access_raster_school <- st_rasterize(access_poly_school["accessibility"], 
                              nx = 1000, ny = 1000) 
plot(access_raster_school)

#Aeshna’s Cycling to School Accessibility Score Map, 4A

# ggplot(Detroit_streets) +
#   geom_stars(data = access_raster_school) +
#   geom_sf(color = "white", alpha = 0.05) +
#   scale_fill_viridis_c(na.value = NA, 
#                        option="A",
#                        name = "Accessibility score \nof bicycle access to\nschools") +
#   theme_void()

plot(image_read("images/Aeshna_4A.png"))

access_points_school <- st_as_sf(access_raster_school, as_points = TRUE)

ggplot(Detroit_streets) +
  geom_sf(data = access_points_school, aes(color = accessibility), size = 0.1) +
  scale_color_viridis_c(na.value = NA, 
                       option="A",
                       name = "Bicycle access to\nschools") +
  theme_void()

access_poly_school2 <- st_as_sf(access_raster_school, as_points = FALSE, merge = TRUE)

ggplot(Detroit_streets) +
  geom_sf(data = access_poly_school2, aes(fill = accessibility), color = 'gray') +
  scale_fill_viridis_c(na.value = NA, 
                       option="A",
                       name = "Bicycle access to\nschools") +
  theme_void()

access_contours_school <- st_contour(access_raster_school, contour_lines = TRUE, 
                              breaks = c(0,1,2,3,4,5,6,7,8))

ggplot(Detroit_streets) +
  geom_sf(color = "gray", alpha = 0.2) +
  geom_sf(data = access_contours_school, aes(color = accessibility), fill = NA) +
  scale_color_viridis_c(na.value = NA, 
                       option="A",
                       breaks = c(0,1,2,3,4,5,6,7,8),
                       name = "Bicycle access to\nschools") +
  theme_void()

#Detroit Libraries

#Keana’s Libraries Map 0B

# Detroit_library <- oe_read("networks/bbbike_detroit.pbf", 
#                    provider = "openstreetmap_fr", 
#                    download_directory = "networks", 
#                    layer = "points", 
#                    quiet = TRUE) %>%
#   filter(str_detect(other_tags, '"amenity"=>"library')) %>%
#   st_filter(Detroit_city_limits) %>%
#   rename(id = osm_id)
# 
# ggplot(Detroit_streets) +
#   geom_sf(color = 'azure2', alpha = 0.2) +
#   geom_sf(data = Detroit_library, color = "darkorchid1") +
#   coord_sf(crs = michigan_state_plane)  +
#   theme_void()

plot(image_read("images/Keana_0B.png"))

grid <- st_sf(st_make_grid(Detroit_city_limits, 
                           square = FALSE, 
                           n = c(100,100),
                           what = "polygons")) %>%
  st_filter(Detroit_city_limits) 

colnames(grid) <- "geometry"
st_geometry(grid) <- "geometry"

grid <- grid %>%
  mutate(id = seq(1, length(grid$geometry), by=1))

grid_points <- st_centroid(grid)

ggplot(grid) +
  geom_sf() +
  geom_sf(data = Detroit_library, color = "chartreuse") +
  geom_sf(data = Detroit_streets, alpha = 0.2) +
  coord_sf(crs = michigan_state_plane) +
  theme_map()

ttm <- travel_time_matrix(r5r_core = r5r_core,
                          origins = Detroit_library,
                          destinations = grid_points,
                          mode = c("WALK"),
                          departure_datetime = as.POSIXct("15-11-2021 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S"),
                          max_walk_dist = 1609.34,
                          max_trip_duration = 480,
                          verbose = FALSE)

ttm_library <- travel_time_matrix(r5r_core = r5r_core,
                          origins = st_transform(x = Detroit_library, crs = "WGS84"),
                          destinations = st_transform(x = grid_points, crs = "WGS84"),
                          mode = "BICYCLE",
                          departure_datetime = as.POSIXct("15-11-2021 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S"),
                          max_walk_dist = 1000,
                          max_trip_duration = 480,
                          verbose = FALSE)

tt_wide_library <- ttm_library %>%
  pivot_wider(names_from = fromId, 
              names_prefix = "from", values_from = travel_time) %>%
  rename(id = toId) %>% 
  merge(grid) %>%
  replace(is.na(.), 999) %>%
  rowwise() %>%
  mutate(from_any = min(c_across(starts_with("from")), na.rm = TRUE))

st_geometry(tt_wide_library) <- "geometry"

tt_wide <- ttm %>%
  pivot_wider(names_from = fromId, 
              names_prefix = "from", values_from = travel_time) %>%
  rename(id = toId) %>% 
  merge(grid) %>%
  replace(is.na(.), 999) %>%
  rowwise() %>%
  mutate(from_any = min(c_across(starts_with("from")), na.rm = TRUE))

st_geometry(tt_wide) <- "geometry"

#Keana’s Cycling to Library Map, 2B

# ggplot(Detroit_streets) +
#   geom_sf(data = tt_wide_library, 
#           aes(fill = from_any), 
#           color = NA) +
#   geom_sf(alpha = 0.1) +
#   scale_fill_gradient2(low = "green", mid = "yellow", high = "red", 
#                        midpoint = 30,
#         name = "bike\ntime to the\nnearest library\n(minutes)") +
#   coord_sf(crs = michigan_state_plane) +
#   theme_map()

plot(image_read("images/Keana_2B.png"))

#Keana’s Library Isochrone, 1B

# iso_pallete <- wes_palette("Darjeeling1", n = 5)
# 
# iso10min <- tt_wide_library[tt_wide_library$from_any < 11,] %>%
#   st_union()
# 
# iso20min <- tt_wide_library[tt_wide_library$from_any < 21,] %>%
#   st_union()
# 
# iso30min <- tt_wide_library[tt_wide_library$from_any < 31,] %>%
#   st_union()
# 
# ggplot(Detroit_streets) +
#   geom_sf(data = iso30min, 
#           aes(fill = "Area within 30 minutes"), 
#           color = NA) +
#   geom_sf(data = iso20min, 
#           aes(fill = "Area within 20 minutes"), 
#           color = NA) +
#   geom_sf(data = iso10min, 
#           aes(fill = "Area within 10 minutes"), 
#           color = NA) +
#   geom_sf(alpha = 0.1) +
#   scale_fill_manual(values = c(iso_pallete[1], 
#                                iso_pallete[3],
#                                iso_pallete[5]),
#         name = "bike\ntime to the\nnearest library\n(minutes)") +
#   coord_sf(crs = michigan_state_plane) +
#   theme_map()

plot(image_read("images/Keana_1B.png"))

library_grid <- grid %>%
  mutate(num_library = lengths(st_covers(grid, Detroit_library)))

library_points <- st_centroid(library_grid)

#Alex’s Libraries Reachable by Cycling Map, 3B

# Detroit_library_access <- accessibility(r5r_core,
#                         origins = library_points,
#                         destinations = library_points,
#                         mode = "BICYCLE",
#                         opportunities_colname = "num_library",
#                         decay_function = "step",
#                         cutoffs = 11,
#                         departure_datetime = as.POSIXct("15-11-2021 14:00:00",
#                                  format = "%d-%m-%Y %H:%M:%S"),
#                         max_walk_dist = 5000,
#                         time_window = 120,
#                         percentiles = 50,
#                         verbose = FALSE) %>%
#   mutate(id = as.numeric(from_id)) %>%
#   merge(grid)
# 
# st_geometry(Detroit_library_access) <- "geometry"
# 
# ggplot(Detroit_library_access) +
#   geom_sf(aes(fill = accessibility), color = NA) +
#   scale_fill_viridis_c(name = "Libraries\nwithin a 20-minute\nbicycle ride") +
#   coord_sf(crs = michigan_state_plane) +
#   theme_void()

plot(image_read("images/Alex_3B.png"))

Detroit_library_access <- accessibility(r5r_core,
                        origins = library_points,
                        destinations = library_points,
                        mode = "BICYCLE",
                        opportunities_colname = "num_library",
                        decay_function = "exponential",
                        cutoffs = 11,
                        departure_datetime = as.POSIXct("15-11-2021 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S"),
                        max_walk_dist = 5000,
                        time_window = 120,
                        percentiles = 50,
                        verbose = FALSE) %>%
  mutate(id = as.numeric(from_id)) %>%
  merge(grid)

st_geometry(Detroit_library_access) <- "geometry"

ggplot(Detroit_library_access) +
  geom_sf(aes(fill = accessibility), color = NA) +
  scale_fill_viridis_c(name = "Accessiblity score") +
  coord_sf(crs = michigan_state_plane) +
  theme_void()

st_write(Detroit_library_access, 'Detroit_library_access.geojson', append=FALSE, quiet=TRUE )

access_poly_library <- st_read("Detroit_library_access.geojson", quiet=TRUE)

access_raster_library <- st_rasterize(access_poly_library["accessibility"], 
                              nx = 1000, ny = 1000) 

Alex’s Cycling to Library Accessibility Score Map, 4B

# ggplot(Detroit_streets) +
#   geom_stars(data = access_raster_library) +
#   geom_sf(color = "white", alpha = 0.05) +
#   scale_fill_viridis_c(na.value = NA, 
#                        option="A",
#                        name = "Accessibility score \nof bicycle access to\nlibraries") +
#   theme_void()

plot(image_read("images/Alex_4B.png"))

stop_r5(r5r_core)